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ABSTRACT 

Understanding  turbulence  degradation  of  electromagnetic 
wave  propagation  is  essential  for  efficient  operation  of  laser 
weapons,  target  designators,  and  imaging  systems.  Random 
atmospheric  refractive  index  inhomogeneities  alter  the  phase 
and  amplitude  of  electromagnetic  waves. 

This  thesis  attempts  to  model  atmospheric  turbulence 
effects  by  using  filtered  Gaussian  phase  screens  to  represent 
the  random  nature  of  refractive  index  changes.  The  simulation 
uses  two-dimensional  512  x  512  fast  Fourier  transform  (FFT) 
techniques  with  extended  Huygens-Fresnel  principles  performed 
on  a  desk  top  computer. 

Simulation  verification  was  accomplished  by  comparing 
calculated  and  theoretical  spatial  coherence  lengths,  ^>  o  . 
Phase  only  screens  produced  coherence  lengths  that  were  30% 
larger  than  theoretical  values.  By  using  random  phase  and 
amplitude  screens,  the  calculated  coherence  lengths  agreed  to 
within  3%  of  theoretical  values.  Saturation  of  the  normalized 
intensity  variance,  a2/!2,  occurred  for  increasing  turbulence 
using  a  single  phase-amplitude  screen. 
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I .   INTRODUCTION 

Turbulence  in  the  atmosphere  has  a  detrimental  effect  on 
the  propagation  of  electromagnetic  waves  at  optical 
wavelengths.  Minute  changes  in  the  refractive  index  of  the 
atmosphere  along  the  path  of  the  wave  cause  perturbations  in 
the  wavefront  in  both  intensity  and  phase.  Consequently, 
random  atmospheric  density  fluctuations  degrade  the 
performance  of  optical  imaging  systems.  Understanding  this 
phenomenon  is  essential  for  future  endeavors  in  laser  weapons 
and  designators,  as  well  as  imaging  systems. 

This  thesis  attempts  to  model  the  effects  of  atmospheric 
turbulence  using  the  extended  Huygens-Fresnel  principle.  This 
theory  considered  both  Fraunhofer  and  Fresnel  propagation 
although  the  simulation  focused  on  a  single  phase  screen  using 
the  Fraunhofer  approximation. 

A  Rytov  approximation  filtered,  Gaussian  distributed, 
random  phase  and  amplitude  screen  simulated  the  effects  of 
turbulence  on  an  optical  beam.  The  magnitude  of  turbulence 
of  this  screen  was  proportional  to  the  index  of  refraction 
structure  parameter,  Cn2. 

The  simulation  code  used  two-dimensional  fast  Fourier 
Transform  (FFT)  algorithms  extensively  to  model  optical  wave 
propagation  in  the  Fraunhofer  regime. 


Accuracy  of  the  simulation  was  checked  by  comparison  of 
calculated  and  theoretical  coherence  lengths,  Do  ,  which  were 
obtained  from  the  e_1  point  of  the  atmospheric  mutual 
coherence  function  (MCF) .  Another  check  verified  that 
intensity  variance  saturation  occured  for  an  increase  in 
turbulence,  as  measured  experimentally  and  predicted  by 
theory. 


II .   BACKGROUND 

A.  STATISTICAL  DESCRIPTION  OF  ATMOSPHERIC  TURBULENCE 
1.   Stationarity ,  Homogeneity,  and  Isotropy 

Most  theorems  involving  random  processes  require  that 
the  functions  satisfy  certain  restrictions.  "Stationarity" 
implies  that  the  mean  value  of  the  random  function  does  not 
change  with  time.  "Homogeneity"  implies  that  translations  in 
x,  y,  and  z  coordinates  (a  Galilean  coordinate  transformation) 
do  not  affect  variables  of  the  random  function.  "Isotropy" 
means  that  the  random  function  is  independent  of  any 
rotational  coordinate  transformation.   [Ref .  1] 

The  last  two  conditions  imply  that  the  statistical 
representation  of  the  random  function  for  two  points  f (ri)  and 
f  (F2)  depend  only  on  the  magnitude  of  their  separation 

I  ri  -  r2  I .   [Ref.  1] 

If  these  three  restrictions  only  hold  for  a  local 
volume  of  space,  L3  ,  the  random  field  is  said  to  be  locally 
stationary,  locally  homogeneous,  and  locally  isotropic  for 

|  Fi  -  Fa  |  <  L. 

Turbulence  of  the  atmosphere  is  a  process  that  can  be 
best  described  by  a  random  function.  Minimal  violations  of 
stationarity,  homogeneity,  and  isotropy  occur  if  the 
neighborhood  is  sufficiently  small,  less  than  some  outer  scale 
length  L.   [Ref.  1] 
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2.  The  Structure  Function 

Because  of  their  nature,  random  functions,  f ,  can  best 
be  described  stochastically.  The  simplest  moment  of  f  is  the 
mean  value,  <f>. 

Another  common  parameter  used  to  represent  a  random 
function  is  the  variance,  o2 ,  defined  as 

o2  =  <  (  f(r)  -  <f(r)>  )2  >.  (2.1) 

In  a  physical  random  field,  such  as  atmospheric 
turbulence,  the  mean  value,  <f (r)>,  is  difficult  to  ascertain 
because  of  trends  in  the  mean  as  well  as  variations  in  spatial 
position,  vertically  and  horizontally.  Therefore,  the  mean 
and  variance  are  awkward  parameters  for  describing  a  real, 
random,  atmospheric  field.  To  resolve  this  difficulty, 
Kolmogorov  introduced  the  structure  function 

Df(r)  =  Df(r2-Fi)  =  <[f(Fi)  -  f(r2)]2>,  (2.2) 

which  resembles  the  variance  at  first  glance.  [Ref.  2]  The 
structure  function  is  an  ensemble  average  taken  over  all 
possible  point  pairs  ri  and  F2 .   [Ref.  2] 

Through  dimensional  analysis,  and  assuming  isotropy, 
Kolmogorov  deduced  a  relation  between  the  structure  function 
and  the  distance  between  the  sampling  points, 

Df  (r)  =  Cf2  r2'3  (2.3) 


which  holds  over  a  limited  volume,  called  the  inertial 
subrange.  [Ref .  2]  Although  developed  for  the  velocity  field 
v,  Equation  2.3  also  holds  for  certain  atmospheric  functions 
called  conservative  passive  additives. 

Conservative  passive  additives  are  quantities  that 
play  no  part  in  the  turbulence  of  the  medium.   Temperature  is 
not  such  a  quantity,  but  potential  temperature,  which  corrects 
for  the  adiabatic  change  in  temperature  with  altitude,  is. 
Therefore,  the  relation 

Dt  (r)  -  Ct2  r2'3  (2.4) 

holds  within  the  inertial  subrange  where  T  is  the  potential 
temperature . 

In  examining  problems  of  optical  propagation,  the 
index  of  refraction  is  the  essential  parameter.  To  express 
the  index  of  refraction  as  a  conservative  passive  additive, 
it  is  useful  to  manipulate  the  following  relation 

n  =  1  +  77.6(1  +  7.52  x  10"3/A2)  (P/T)  x  10"8    (2.5) 

where  X  is  the  wavelsngth  of  the  radiation  in  pm,  P  is  the 
atmospheric  pressure  in  mb,  and  T  is  the  atmospheric 
temperature  in  K.  [Ref.  3]  The  differential  of  Equation  2.5 
is 

dn  =  -  79  x  10-*  P  dT/T2  (2.6) 


for  red  light  (X  =  0.6  um)  after  ignoring  the  dP  term 
(isobaric  turbulence) .  Since  potential  temperature  is  a 
passive  additive,  the  index  of  refraction  is  also.  Using 
Equation  2.6  to  express  the  index  of  refraction  structure 
parameter  in  terms  of  the  temperature  structure  parameter 
gives 

D„(r)  =  C»2  r2'3  =  (79  P/T2  x  10-*)2  d2  r2'3     (2.7) 
[Ref.  3]. 

3.  The  Correlation  Function  and  Power  Spectral  Density 
Expanding  Equation  2.2  gives 

Df(r)  =  <f2(rx)>  -  2<f (ri)f (r2)>  +  <f2(r2)>.     (2.8) 

But  if  we  assume  an  homogeneous  medium,  then  <  f2 (r)  >  is  the 

same  for  any  point  r,  assuming  <f(r)>  =  0.  The  structure 
function  becomes 

Df (r)  =  2  Bf (0)  -  2  Bf (r),  (2.9) 

where  Tatarski  defined  the  correlation  function,  Bf ,  as 

Bf  (r)  =  <  f (ri)  f* (r2)  >  (2.10) 

in  a  stationary,  homogeneous,  and  isotropic  medium.  [Ref.  2] 
Using  stochastic  Fourier-Stielt jes  integrals,  Tatarski 
developed  a  relation  for  the  three-dimensional  power  spectral 
density  function,  *(f),  as  the  Fourier  transform  of  the 
correlation  function.  Writing  this  in  freqency  ( f)  form  as 
opposed   to  radian  freqency  (K)      form  preferred  by  Tatarski, 


exp  (-2nifT)  B(r)  d3r,         (2.11) 


and  conversely, 

Bf(r)  = 


exp  (-2nif--f)  *(f")  d3  f.  (2.12) 


[Ref.  2] 

Using  the  fact  that  the  correlation  function  is  even,  and 
using  relations  (2.7)  and  (2.9),  the  structure  function  for 
the  Kolmogorov  power  spectral  density  is 

*B(f)  =  1.303  Cn2  f-ii/a  ,  (2.13) 

for  the  Kolmogorov  power  spectral  density.   This  equation  is 
analogous  to 

4>„(/D  =  0.033  Cn2  Jr-1*'3,  (2.14) 

an  expression  developed  by  Tatarski  in  K   space  that  is  seen 
in  many  publications.   [Ref.  2] 

B.   ELECTROMAGNETIC  PROPAGATION  THROUGH  TURBULENCE 
1 .   The  Wave  Equation 

Central  to  any  study  of  electromagnetic  propagation 
are  the  four  equations  of  Maxwell,  presented  here  in  Gaussian 
units 

V  •  H  =  0,  (2.15) 

V  x  H  =  -ikn2E,  (2.16) 


V  .  (n2E)  =  0,  (2.17) 

V  x  E  =  ikfi,  (2.18) 

for  the  assumptions  that  the  medium  has  zero  conductivity, 
unit  magnetic  permeability,  and  sinusoidal  dependence. 
Taking  the  curl  of  Equation  2.18  gives 

V  x  (VxD  -  V  x  (ikH).  (2.19) 

Using  a  vector  identity  and  Equation  2.16,  Equation  2.19 
becomes 

-V2E  +  V(V  •  E)  =  k2n2E.  (2.20) 

Expanding  Equation  2.17  gives  a  relation  for 

V  •  E,  which  when  inserted  into  Equation  2.20  gives  the  wave 

equation 

V2E  +  k2n2E  +  2V[E  •  V(log  n) ]  =  0.        (2.21) 

If  the  wavelength  of  the  electromagnetic  wave  is  small 
compared  to  the  dimensions  of  the  refractive  inhomogeneities , 
then  the  2V[E  •  V(log  n) ]  term  is  negligible.  In  this  case, 
the  wave  equation  reduces  to  the  simple  form 

V2E  +  k2n2E  =  0.  (2.22) 

[Ref.  3] 

2.  The  Born  Approximation 

One  method  of  solving  the  wave  equation,  employed  by 
both  Tatarski  [Ref.  2]  and  Clifford  [Ref.  3],  is  the  method 
of  small  perturbations  or  the  Born  approximation.  This  method 
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expands  the  electric  field  and  index  of  refraction  as  a  power 
series 

E  =  Eo  +  Ei  +  •  •  •  ,  (2.23) 

n  ■  1  +  xu  +  •••,  (2.24) 

where  Eo  is  a  constant  in  time,  and  Ei  and  m  are  time 
varying. 

Inserting  these  two  relations  into  Equation  2.22  and 
equating  same  order  terms  gives 

V2Eo  +  k2Eo  =  0,  (2.25) 

V2Ei  +  k2Ei  +  2mk2Eo  =  0,  (2.26) 

ignoring  all  second  and  higher  order  terms. 

Assuming  that  the  electromagnetic  wave  propagates  in 
the  z-direction  and  E  o  =  exp(ikz)  ,   then  Equation  2.26  becomes 

V2Ei  +  k2Ei  +  2mk2exp(ikz)  =  0.         (2.27) 

Clifford  solved  this  wave  equation  for  Ei  by  utilizing 
a  Green's  function.  The  solution  is  the  convolution  of  a 
plane  wave  Green's  function  with  the  source  term,  shown  here 
in  scalar  form 

1   C      exp[ik|?-r'  |] 
Ei(r)  =  —   /  [2k2  ni(r')  exp(ikz')]  d3r, 

4n  Jv  |r-r' | 

(2.28) 

where  r  is  the  position  of  a  specific  point  in  the  image  plane 
and  r'  is  the  postion  of  specific  source  point.    [Ref .  3] 


Figure  (2.1)  shows  the  coordinate  system  used  to  express  the 
solution  to  the  wave  equation  for  propagation  between  planes 
z'  and  5'    to  z  and  p  . 

For  normal  propagation  situations,  jz-z*|  =  R  >>  \  f>  ~ 
p%\    so  that  Equation  2.28  expands  into  the  form 

1  r  exptikR  [1  +  \p-  5*  \*/2R*    +  •••]) 

Ei(r")  =  —  /  

4njv      R  [1  +  |^-5'|2/2R2  +  •••] 

X   [2k2  m(r')  exp(ikz')]  d3r,        (2.29) 

using  the  binomial  expansion.  Dropping  those  terms  consistent 
with  the  Fresnel  approximation  in  Huygens-Fresnel  optics 
theory,  this  equation  reduces  to 


(2.30) 

which  is  known  as  the  Fresnel  diffraction  formula.   [Ref.  3] 
3 .  The  Rytov  Approximation 

Another  method  of  solving  the  wave  equation,  also  used 
by  both  Tatarski  and  Clifford,  is  the  method  of  smooth 
perturbations  or  the  Rytov  approximation.  This  technique 
assumes  that  the  electric  field  is  of  the  form 

E  =  exp(^)  =  exp(X+iS)  =  A  exp(iS).       (2.31) 
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Figure  2.1   Aperture  and  Image  Plane  Geometry 
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Using  power  series  expansions  for  the  electric  field, 
Equations  2.25  and  2.26  still  hold.  Placing  Equation  2.31 
into  these  two  gives 

(V*E)/E  +  k*n*(r)  -  0.  (2.32) 

Tatarski  shows  that  this  method  gives  the  sane  results  as  the 
method  of  small  perturbations  but  the  range  of  validity  is 
larger  than  for  the  Born  approximation.  Equation  2.26  holds 
as  long  as  V*E  is  small  compared  to  A.   [Ref .  2] 

C.   HUYGBMS-FRESHEL  THEORY 

Another  way  of  obtaining  a  relation  for  the  electric  field 
amplitude  at  the  observation  plane,  is  by  extending  the  ideas 
of  the  Huygens-Fresnel  theory  to  include  a  random  phase  term. 
The  Huygens-Fresnel  theory  offers  this  solution  for  the 
electric  field  seen  in  the  observation  plane 

E(r)  =  -  ik  I   E(r')   exp(ik|r-r' |)  dp  '.       (2.33) 
2n  J s  |r-f '  | 

[Ref.  5] 

Lutomirski  and  Yura  develop  an  extension  of  the  Huygens- 
Fresnel  idea  which  incorporates  a  random  phase  term  which  is 
equivalent  to  the  form  used  in  the  Rytov  approximation, 
exp  (  fy  )  ,  where  f  is  a  complex  phase.  [Ref.  4]  This  extended 
integral  is 

E(r)  =  -  ik   /   E(r')  exp(^(f))  exp(ik|f-f '  |)  dp   '. 

2n  JB  !?~?,| 

(2.34) 
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In  the  geometrical  optics  limit,  j   (r)  becomes  the  optical 
phase 

V  (r)  =  ik  /  m(z)  dz.  (2.35) 

Expanding  ty    into  a  power  series,  Equation  2.35  reduces  to 
Equation  2.30  as  developed  by  Tatarski  and  Clifford. 

D.   MUTUAL  COHERENCE  FUNCTION 

In  trying  to  get  a  quantitative  measure  for  turbulence, 
Lutomirski  and  Yura  [Ref.  4]  examined  the  average  intensity 
of  the  wave  at  the  observation  point.  The  average  intensity 
at  r  is  defined  as 

<I(r)>  =  <  E(ri)  EMr2)  >.         (2.36) 

From  Equation  2.34: 

<I(r)>  =  <  (k*/4n2)  //  E(Fi  •  )E*  (r* ■ )  exp[^(fi)  +  ^*(r«)] 

X      exp[ik(  jri-ri  '  \    -    |r2-r2'|)3    d  p  i      dp  2      > 

|ri-ri'|     Ifi-rVI 

(2.37) 

Picking  out  only  the  time-dependent  portion  of  this  result, 

the  entire  integral  reduces  to  evaluating 

<  exp[^(ri)  +  ^*(r2)]  >.  (2.38) 

This  term  is  the  atmospheric  mutual  coherence  function  (MCF) 

which  contains   the  elements  of  the  propagation  through 

turbulence . 

Lutomirski  and  Yura  [Ref.  4]  relate  the  atmospheric  MCF 

to   the  wave   structure   function,    D,     by   the  relation 
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MCF(^  )  =  <exp[^(ri)  +  ^Mr2)]>  =  exp[-D(^  )/2]     (2.39) 
where  Dip  )    =   Dx  ( p  )    +  Ds  ( p  )    from  the  Rytov  relation  <^  -   X+iS 


or 


MCF 


(p)    =  exp[-(/7/^o)s/3]  ,  (2.40) 


where 


p0    =  [  1.46  k2   /"  Cn2  dz  ]-3'9.        (2.41) 


0 


and  k  is  the  wavenumber .  Cn2  is  the  refractive  index 
structure  parameter  defined  by  Equation  2.7,  and  R  is  the 
distance  of  propagation  through  the  atmosphere.  The  coherence 
length,  Do,  is  defined  as  the  distance  where  transverse 
spatial  coherence  of  the  wave  drops  to  the  e-1  .   [Ref.  1] 

E.   THE  DIFFRACTION  INTEGRAL 

Like  Lutomirski  and  Yura,  Roberts  [Ref.  6]  begins  with  the 
Huygens-Fresnel  integral  for  the  propagation  of  light  waves 

after  making  the  Fresnel  approximation  (jz-z'j  =  R  >>  Jo-5' | ) 


E(r)  =  -ik    /  dp    '    E(r')  exp{  ik  [R2+ |  Z)  '-p   |2]1/2.   (2.42) 
2nR 


In  this  expression,  k  is  the  wavenumber  and  the  notation  is 
the  same  as  shown  in  Figure  (2.1). 

Employing  the  standard  practice  of  expanding  by  the 

binomial  expansion  and  neglecting  smaller  terms  gives 

E(r)  =  -ik   exp[-ikR]  /  dp   '  E(r')  exp{ik|  p'-p   |2  I  • 

2nR  Is  2R 

(2.43) 
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Since  the  exponential  term  outside  the  integral  doesn't  affect 
intensity  measurements,  it  can  be  neglected.  Expanding  the 
quadratic  term  gives  the  relation 

E(r)  =  -ik_   /  dp    '  E(r')  expfik  [  p  '2-2  p  '  •  p  +    p2]\ 
2nR   /s  '  2R    r     il' 

=  -ik_  exp[ik  p'2]   dp'    E(r')  exp[ik 2-]  exp[-ik>5  ■  5  '1  . 

2nR      2Rr     /  s  2R         ~^ 

(2.44) 

This  is  called  the  convolution  form  of  the  Fresnel  diffraction 

integral  which  differs  from  the  Fraunhofer  approximation  only 

by  the  quadratic  phase  term  exp[ik/j_2].   [Ref.  6] 

If  we  denote  the  Fourier  transform  of  E(r)  as  E{f),    then 


***       •- 


E(f)    =     /  d^   exp(-2nif-/5  )  E(r).  (2.45) 


Inserting  E(r)  as  defined  in  Equation  2.43,  gives 


E(f)    =  -ik    dp   '    E(r')    dp     exp(-2nif-p   )    exp{  ik  |  p   '  -  p   |2  . 

2nR  Is  f a  2R 

(2.46) 


Changing  variables  to  q  =  p   '  -  D     ,     the  electric  field 
spectrum  becomes 

E(f)    =  -ik    /  dp'    E(r')  exp(-2nif-^  ) 

2nR  Js  ' 

X    /dq  exp(-2nif-q)  exp(ik  q2).        (2.47) 

Js  2R 

The  q  integral  is  shown  to  be  the  Fourier  transform  of  the 
Gaussian  function 
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2niR  exp(-2n2iRf2)  .  (2.48) 

k         k 


[Ref .  6]   Substituting  this  function  into  Equation  2.47  gives 


E(f)    =  exp  ( -2niiR  f2  )   /  dp   ■  E(r')  exp(-2nif  >p  )  ,   (2.49) 


which  is  the  inverse  transform  of  what  we  sought, 


E(r)  =   /  df   exp(2nif«r)  exp  (-2n2-iRf2  ) 

k 


X  J  dp    '    E(r')  exp(-2nif-y0  )  (2.50) 

This  is  the  transfer  form  of  the  Fresnel  diffraction  integral. 
[Ref.  6] 

Equations  2.44  and  2.50  are  two  different  forms  of  the 
diffraction  integral  that  are  useful  for  two  types  of  wave 
propagation.  The  convolution  form,  Equation  2.44,  is  best 
suited  for  long  distance  propagation  (Fraunhofer  regime)  since 
R  is  in  the  denominator  of  the  exponential  term.  The  transfer 
form,  Equation  2.50,  is  best  suited  for  short  distance 
propagation  since  R  is  present  in  the  numerator  of  the 
exponential  term.   [Ref.  6] 
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III.  NUMERICAL  SIMULATION 

A.   SIMULATION  HARDWARE  AND  SOFTWARE  SPECIFICS 

This  numerical  simulation  modeled  the  propagation  of  an 
electromagnetic  wave  through  a  random  turbulent  medium.  The 
procedures  for  this  simulation  required  the  creation  and 
manipulation  of  multiple  two-dimensional,  512  by  512,  single 
precision,  phase  screens,  requiring  a  minimum  memory  of 
several  megabytes. 

A  COMPAQ  Deskpro  80386-20  computer  configured  with  a  16 
megabyte  RAM  and  64  megabyte  hard  drive  met  hardware 
requirements.  The  computer  also  had  EGA-VGA  graphics  and  both 
HP  Laser  Jet  II  and  Panasonic  KX-P1092i  printers.  A  Weitek 
1167  math  coprocessor  increased  execution  speeds  of  the  20  Mhz 
80387  coprocessor  by  a  factor  of  3-4. 

The  software  requirements  for  writing  the  simulation  code 
were  met  with  the  MicroWay  32  bit  NDP  FORTRAN-386  compiler 
and  the  Phar  Lap  operating  system  extender.  It  was  chosen 
over  the  Science  Applications  International  Corporation's  SVS 
Fortran-386  compiler  which  had  10-20%  faster  program  execution 
times  but  numerous  bugs  present  in  the  graphics  routines. 

Unfortunately,  due  to  limitations  of  version  1.4e  of  the 
NDP  Fortran-386  compiler,  the  full  graphic  capabilities  of  VGA 
were  not  supported,   so  EGA  graphics  were  used  for  screen 
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displays.   Hardcopy  graphs  and  plots  were  accomplished  using 
Grapher  and  Surfer  programs  from  Golden  Software. 

B.   REQUIRED  ALGORITHMS 

For  the  successful  operation  of  such  a  model  both  the 
random  number  generator  and  the  two-dimensional  fast  Fourier 
transform  (FFT)  routines  needed  speed  and  accuracy.   These 
will  be  discussed  in  detail  in  the  following  sections. 
1 .  Random  number  generator 

The  effects  of  turbulence  on  the  propagation  of  an 
electromagnetic  wave  are  of  a  random  nature.  To  emulate  this 
random  nature  required  a  pseudo-random  number  generator 
algorithm.  The  minimum  criteria  for  a  quality  random  number 
generator  was  that  the  products  be:  (1)  sufficiently  random, 
(2)  reproducible,  and  (3)  rapidly  obtained. 

The  NDP  FORTRAN  compiler  did  not  come  equipped  with  a 
random  number  generator,  so  an  algorithm  of  the  linear 
congruent  type  especially  designed  for  32  bit  machines  was 
implemented  as  a  subroutine  attached  to  the  main  program. 
This  had  the  benefit  of  not  binding  this  portion  of  the  code 
to  a  specific  machine  setup.  This  random  number  generator 
algorithm  came  from  Dr.  Milne,  who  obtained  it  from  class 
notes  prepared  by  Dr.  Harrison  [Ref.  7].  This  algorithm 
appears  as  subroutine  RAN   (Appendix  A)  in  the  simulation  code. 

The  first  criteria,  that  the  generator  produce  numbers 
sufficiently  random,  was  hard  to  quantify.   A  series  of  tests 
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exercised  the  statistical  independence  of  succeeding  numbers. 
There  are  many  such  tests  that  exist,  the  most  common  of  which 
is  the  MChi-Square"  test.  The  following  is  a  list  of  all 
tests  used  to  check  the  quality  of  the  random  number  generator 
algorithm: 

1.  Frequency  Test.  This  tests  the  uniformity  of  the 
sequence  of  numbers. 

2.  Serial  Test.  This  tests  the  two-dimensional  uniformity 
of  the  sequence  of  numbers . 

3.  Runs  Up  and  Down  Test.  This  tests  to  see  how  long 
sequences  of  random  numbers  are  that  either  go 
successively  up  or  down. 

4.  Laqqed-Product  Test.  This  tests  for  correlations 
between  random  numbers  Ri  and  Ri  +  j  where  j  is  the  given 
lag  value.   A  lag  of  3  is  especially  difficult  to  pass. 

5.  Repeat  Test.  This  simple  qualitative  check  is  to 
determine  how  long  the  sequence  of  random  numbers 
becomes  before  the  first  number  is  repeated. 

6.  Scatter  Plot  Test.  This  is  another  qualitative  test 
where  random  numbers  plotted  on  a  x-y  plane  are  visually 
checked  for  correlation. 

These  six  tests  were  performed  for  a  variety  of  seed 
values  to  determine  which  gave  the  best  performance.  Results 
of  the  first  five  tests  for  two  good  seeds  can  be  found  in 
Table  3.1. 

Choosing  the  same  initial  seed  value  gave 
reproducibility  of  the  random  numbers.  The  algorithm's  simple 
two  lines  of  code  assured  rapid  execution. 
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TABLE  3.1 


RANDOM  NUMBER  GENERATOR  TEST  DATA 


This  table  provides  the  results  of  statistical  tests  of 
the  random  number  generator  algorithm.  The  X2  values  came 
from  the  Chi-Square  table  in  Bevington  [Ref.  11] 

The  results  of  the  Lagged  Product  test  correspond  to  the 
theoretical  mean,  ut  ,  the  calculated  mean,  p,  the  theoretical 
standard  deviation,  ot  ,  and  the  calculated  standard  deviation. 


Test 

Seed 

Degrees  of 
Freedom 

X* 

Probability 

Frequency 

45739 
94377 

9 
9 

4.0 
5.1 

91  % 
82  % 

Serial 

45739 
94377 

20 
20 

16.6 
10.2 

70  % 
96  % 

Runs 

Up  and 

Down 

45739 
94377 

6 
5 

11.7 

2.4 

2  % 
80  % 

Test 

Seed 

ut 

P 

Ot 

o 

Lagged 
Product 

45739 
94377 

0.250 
0.250 

0.230 
0.265 

0.088 
0.088 

0.092 
0.086 

2.   Fast  Fourier  Transform 

Since  the  most  often  used  algorithm  in  the  simulation 
was  the  Fast  Fourier  Transform  (FFT) ,  efficiency  was  critical. 
Two  different  routines  were  compared. 

The  first  FFT  routine  considered  was  a  package  offered 
by  NDP  of  one-  and  two-dimensional  machine  language  FFT 
routines  designed  for  easy  linking  by  the  NDP  FORTRAN 
compiler.  This  package  had  routines  available  for  real  and 
complex  arrays. 
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The  second  FFT  routine  considered  was  a  simple  real 
array  subroutine  coded  by  Dr.  Walters  from  a  BASIC 
demonstration  package  provided  by  Infotek. 

Each  FFT  routine  was  compared  for  numerical  output  and 
speed  of  operation  for  the  one-dimensional  case.  As  expected, 
all  routines  were  virtually  identical  in  numerical  output. 
However,  speed  of  operation  differed  widely,  as  seen  in  Table 
3.2. 


TABLE  3.2 


FFT  SPEED  OF  OPERATION  COMPARISON 


This  table  contains  results  of  FFT  execution  time 
comparison  between  NDP  machine  language  routines  and  Dr. 
Walters'  subroutine.  Since  NDP  FORTRAN  does  not  have  a  timing 
function  with  sufficient  accuracy,  execution  times  listed  are 
mean  times  found  from  several  runs,  timed  by  stopwatch. 


Array 
Size 

NDP  Routines 
CFFT           RFFT 

Walters '  Routine 
w/o   &   w/  Weitek 

1-D  field  size 

time (s) 

time (s) 

time (s) 

time (s) 

213 

2»< 

213 

216 
217 

2.60 

5.27 

10.94 

22.78 

47.61 

1.45 

2.88 

5.84 

12.05 

25.02 

2.84 

5.86 

12.27 

25.90 

54.53 

1.20 

2.32 

4.94 

10.26 

21.55 

The  differences  in  run  time  for  the  two  NDP  machine 
language  routines  were  due  to  the  CFFT  routine  being  more 
cumbersome  in  converting  repeatingly  between  real  and  complex 
arrays.  The  compiled  subroutine  was  significantly  faster  than 
both  CFFT  and  RFFT  NDP  machine  language  routines  unless  the 
Weitek  was  deactivated.   Using  a  subroutine  enclosed  within 
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the  main  program  had  the  additional  benefit  of  not  being 
machine  dependent.  Therefore,  it  was  decided  that  Dr. 
Walters'  subroutine  would  be  used  for  all  FFTs  needed  in  the 
simulation  code.  Conversion  of  this  subroutine  for  two- 
dimensional  use  was  relatively  easy.  The  NDP  two-dimensional 
FFT  routines  never  worked  correctly. 

A  FFT  "breaks  down"  a  function  into  its  cosine  and 
sine  constituents  (real  and  imaginary  terms,  respectively). 
The  range  of  spatial  frequencies  represented  in  a  FFT  depend 
on  the  array  size  chosen.  The  minimum  frequency  is  f«m  =  1/W 
and  the  maximum  frequency  is  f«ax  =  n/W,  where  W  is  the 
physical  width  of  the  array  and  n  is  the  number  of  sampling 
points  across  the  array.  When  energy  at  frequencies  near 
these  cutoff  values  appeared,  problems  developed  due  to  the 
discrete  nature  of  the  FFT. 

One  problem  exists  in  representing  spatial  frequencies 
smaller  than  fain  due  to  the  high  amplitude,  low  frequency 
"tilt"  of  the  f  -11'3  filtered  wavefront.  Since  a  tilted 
straight  line  approximates  a  cosine  or  sine  function  for  only 
a  small  region,  it  helps  to  make  the  wavefront  smaller  than 
the  FFT  array  size.  Hence,  to  help  represent  lower 
frequencies  in  the  simulation  coding  accurately,  an  aperture 
array  was  a  user  chosen  variable  smaller  than  the  FFT  array 
size . 

Another  problem  that  manifests  itself  is  "aliasing", 
which  occurs  when  high  frequency  features  of  a  function  are 

22 


missed  by  having  too  few  sample  points.  This  means  that 
important  information  for  the  function  present  in  the 
frequencies  greater  than  faax  are  folded  back  and  appear  at 
lower  frequencies.  Choosing  the  distance  between  sample 
points  smaller  than  the  coherence  length  of  the  propagating 
electromagnetic  wave,  found  theoretically  by  Equation  2.39, 
alleviated  this  problem. 

C.   SIMULATION  DESCRIPTION 

1.   User  defined  quantities 

For  operation  of  the  simulation  program,  the  user  had 
the  option  of  varying  several  quantities  directly  from  the 
keyboard: 

1.  Array  size.  This  can  be  varied  up  to  a  maximum  of  512 
by  512  points. 

2.  Aperture  sampling  points.  This  can  be  varied  up  to  a 
maximum  of  512  by  512  points.  Ideally  it  should  be 
smaller  than  the  array  size  to  avoid  "tilt"  effects,  but 
large  enough  to  avoid  "aliasing"  effects.  A  size  one- 
fourth  to  one-half  the  array  size  works  best. 

3.  Aperture  shape.  A  choice  of  square  or  circular  aperture 
was  provided. 

4.  Seed .  The  pseudo-random  number  generator  required  a 
five  digit  input  seed.  Table  3.1  lists  two  quality 
seeds.   All  data  was  obtained  using  the  seed  94377. 

5.  Turbulence  parameter.  Any  value  of  Cn2  above  zero  can 
be  chosen.   Typical  values  of  Cn2  were  10"17  through 

10-i3  m-2'3. 

6.  Wavelength .  Any  value  of  wavelength  for  the  electro- 
magnetic wave  can  be  chosen.  A  value  of  0.6  urn  was  used 
throughout . 

7.  Distance .  This  was  the  distance  from  the  aperture  to 
the  observer.   This  simulation  used  1  km. 
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8.  Aperture  width.  Adjustable  to  ensure  that  "aliasing" 
effects  do  not  occur  for  a  specific  choice  of  Cn2  . 
Typically,  a  100  mm  aperture  was  chosen  with  a  path 
length  of  1  km. 


For  the  accumulation  of  multiple  data  points  for 
statistical  purposes,  a  slightly  different  version  of  the 
program  read  user  choices  from  a  disk  file. 

2.  Aperture  Mutual  Coherence  Function 

Computing  far-field  diffraction  patterns  of  intensity 
for  an  electromagnetic  wave  through  the  aperture  with  no 
turbulence  gave  the  Mutual  Coherence  Function  (MCF)  of  the 
aperture.  Since  we  assumed  a  point  source  at  infinity,  the 
electric  field  wave  front  was  planar  at  the  aperture.  The 
electric  field  that  passes  through  had  a  value  of  one  within 
the  boundaries  of  the  aperture,  and  a  value  of  zero  outside 
the  boundaries.  The  FFT  subroutine  computed  the  diffraction 
pattern,  the  electric  field  present  at  the  image  plane.  The 
conjugate  square  of  this  electric  field  was  the  intensity 
field  that  would  be  seen  by  an  observer.  Taking  the  inverse 
FFT  of  the  intensity  field  gave  the  aperture  MCF,  which  was 
identical  to  the  autocorrelation  of  the  aperture  function. 

3.  Phase  screen  generation 

Simulation  of  effects  of  turbulence  requires  that  the 
aperture  planar  electric  field  be  multiplied  by  the  phase 
screen.   Two  methods  exist  to  create  the  phase  screens. 
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a.  Classical  method 

In  this  method  n2  Gaussian  distributed  random 
numbers  filled  a  n  by  n  real  array.  Since  the  pseudo-random 
number  generator  RAN  produced  uniformly  distributed  numbers, 
a  subroutine  GAUSS  transformed  them  to  a  Gaussian  distribution 
with  unit  variance,  Knuth  [Ref.  8].  This  subroutine  appears 
in  Appendix  B.  Taking  a  direct  FFT  creates  two  Gaussian 
distributed  arrays,  phaser  and  phasei,  which  are  the  real  and 
imaginary  spectral  amplitudes  of  the  original  random  array. 
Since  these  two  arrays  are  in  spatial  frequency  space, 
multiplication  by  Equation  2.13  occurs.  An  inverse  FFT 
returns  the  filter  to  real  space.  Because  of  the  symmetry  of 
the  filter,  discussed  later,  only  one  array  in  real  space 
results,  containing  all  n2  pieces  of  information  produced  by 
the  subroutine  RAN. 

b.  Martin  and  Flatte  method 

In  this  second  method  suggested  by  Martin  and 
Flatte  [Ref.  9] ,  two  n  by  n  arrays,  called  phaser  and  phasei, 
are  each  filled  with  n2  Gaussian  distributed  random  numbers. 
These  are  considered  the  real  and  imaginary  components  in 
spatial  frequency  space.  Filtering  occurs  as  f_11/3  then  an 
inverse  FFT  returns  the  array  to  real  space.  Because  no 
information  is  lost  in  filtering,  two  arrays  result,  each  with 
n2  pieces  of  information.  Only  one  of  these  arrays  need  be 
used  for  the  following  simulation  code.  This  method  requires 
twice  the  random  numbers  as  the  previous  method  but  GAUSS 
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produces  two  random  numbers  and  half  as  many  two-dimensional 
FFT  steps  are  required.  For  use  in  Fresnel  propagation,  where 
multiple  phase  screens  exist,  this  method  is  certainly  the 
most  efficient. 

4.   Filtering 

Phase  screens  need  to  be  filtered  in  frequency  space 
to  ensure  the  proper  power  law  form.  From  Tatarski  [Ref.  1] 
and  Martin  and  Flatte  [Ref.  9],  the  correct  form  for  the 
filtering  function,  <f>t  ,  is 

*t  =  2n  k2  6z  $„  ,  (3.1) 

where  k  is  the  wavenumber  of  the  electromagnetic  wave,  52  is 
the  slab  thickness,  and  $u  is  the  power  spectral  density 
function.  Using  Equation  2.13  for  the  power  spectral  density 
function  and  the  definition  of  wavenumber,  k  =  2n/A ,  gives 

<t>f    =    1.303     (2n)3     (1A)2    Cn2    52    f  -11'3  (3.2) 

which  is  the  filtering  function  used  in  the  subroutine  FILTER 
in  the  simulation  code.  Arrays  phaser  and  phasei  are  both 
multiplied  by  the  square  root  of  Equation  3.2  to  obtain  a 
factor  for  the  real  and  imaginary  amplitudes  so  that  the 
intensity  has  the  proper  r"-1,/3  power  law. 

Although  this  filtering  scheme  looks  simple,  a  few 
subtle  points  make  the  process  much  more  complicated  than  is 
immediately  evident.  In  fact,  perfecting  this  filtering 
process  constituted  the  major  effort  of  this  thesis. 
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First  of  all,  most  of  the  relations  developed  in  the 
background  references  were  in  terms  of  the  radian  frequency 
K.  The  FFT  routine  used  was  in  terms  of  the  spatial  frequency 
f  =  K/2n.  Consequently,  all  relevant  equations  had  to  be 
expressed  in  this  latter  form,  requiring  the  tracking  down  of 
many  2n  terms. 

Secondly,  due  to  the  two-dimensional  isotropy  of 
turbulence,  the  correct  filtering  scheme  is  of  a  circular 
nature.  This  means  that  the  frequency  seen  in  Equation  3.2 
is  really  f  =  (fx2  +  f"y2)1/2,  which  is  radial  everywhere  but  at 
the  zero  frequency  (DC  term)  where  it  is  has  a  value  of  zero 
[Ref .  9] . 

For  convenience  of  coding,  the  zero  frequency  should 
be  at  the  center  of  the  arrays  phaser  and  phasei.  But  in  the 
FFT  algorithm,  spatial  frequencies  are  arranged  in  a  manner 
that  is  quite  different.  The  zero  frequency  term  is  at  the 
upper  left  corner  of  the  array.  The  Nyquist  "folding" 
frequency  is  present  as  a  cross  that  divides  the  arrays  into 
four  sections.  So,  prior  to  filtering,  a  shuffling  of  data 
points  is  required  as  discussed  by  Brigham  [Ref.  10]  and  as 
seen  in  Figure  3.1.  After  filtering,  an  inverse  shuffle 
returns  all  frequencies  to  their  previous  locations. 

Another  important  consideration  in  the  filtering 
process  was  that  frequency  units  be  correct  according  to  the 
physical  dimensions  of  the  aperture.  This  required  that  the 
frequency  array  have  a  normalization  factor  of   1/NA  where  N 
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Figure  3.1    Graphical  Presentation  Of  Two-dimensional 
FFT  Reorganization  Required  For  Conventional  Viewing. 
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is  the  total  number  of  points  across  the  array,  and  A  is  the 
sampling  interval  (physical  units  between  the  aperture  array 
points) . 

5.   Propagation 

The  Huygens-Fresnel  techniques  formed  the  basis  for 
the  simulation's  treatment  of  optical  propagation.  For 
simplicity  this  simulation  only  considered  Fraunhofer 
propagation.  The  addition  of  subroutines  to  handle  the 
dynamics  of  Fresnel  propagation  should  not  be  too  difficult 
other  than  the  need  for  dynamic  scaling  of  the  numerical  mesh 
to  account  for  diffraction  spreading. 

Martin  and  Flatte  provided  a  procedure  to  simulate  the 
propagation  of  a  wave  in  a  turbulent  medium.   [Ref.  9] 


1.  Add  the  random  turbulence  effects  to  the  electric  field 
distribution  by  meshing  the  aperture  electric  field 
function  with  a  random  phase  term  exp(V): 

E(r')  -  E(r')  exp(^  )  .  (3.3) 

2.  Take  the  Fourier  transform  of  the  result  of  Equation  3.3 
to  obtain 

E{f'  )  =  DFT{E(r'  )  }  .  (3.4) 

3.  Add  a  term  to  the  result  of  Equation  3.4  to  obtain 

E(f)    =    E(f')    expl-aniRf2!  (3.5) 

k 

which  is  the  exact  same  relationship  as  in  Equation 
2.48  that  was  described  by  Roberts  as  the 
convolution  form  of  the  Fresnel  diffraction  integral. 
[Ref.  6] 

4.  Take  the  inverse  Fourier  transform  of  Equation 
3.5  to  get  an  expression  for  the  electric  field 
distribution  at  the  image  plane: 
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E(r)  =  IFTi  E(f)  ]  (3.6) 

If  one  were  taking  a  multiple  screen  approach  to 
solving  for  the  propagation  of  the  wave  as  suggested  by  Martin 
and  Flatte,  the  entire  process  would  be  reiterated  for  as  many 
segments  of  length  R  as  required.  In  this  numerical 
simulation,  only  one  segment  was  used  for  all  propagations. 
Adding  this  iterative  section  to  the  code  would  be  a 
straightforward  modification  if  dynamic  array  sizes  were 
incorporated.   [Ref.  6] 

For  distances  of  propagation  where  the  phase  front 
cannot  be  assumed  planar,  an  additional  quadratic  term  would 
have  to  be  added  at  step  3  for  the  transfer  form  of  the 
Fresnel  diffraction  integral,  as  seen  in  Equation  2.50. 
6.   Atmospheric  Mutual  Coherence  Function 

We  computed  the  atmospheric  mutual  coherence  function 
after  meshing  the  random  phase  screen  with  the  aperture 
electric  field  function  (after  step  1  of  the  Martin  and  Flatte 
method)  .  As  with  the  aperture  MCF ,  the  next  step  was 
computing  the  diffraction  pattern  using  the  FFT  routine,  which 
represents  the  complex  electric  field  that  an  observer  would 
see  because  of  the  aperture  and  the  turbulence.  The  effects 
of  turbulence  tends  to  spread  out  the  diffraction  pattern  so 
that  resolution  would  be  lost  for  any  optical  imaging  system. 
Taking  the  conjugate  square  of  this  diffraction  pattern  gives 
the  intensity  field,  which  negates  the  additional  exponential 
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term  present  in  Equation  3.5.  Taking  the  inverse  FFT  of  the 
intensity  field  gave  the  total  mutual  coherence  function. 
Dividing  the  total  MCF  by  the  aperture  MCF  leaves  the 
atmospheric  MCF. 

7 .   Coherence  Length 

As  stated  previously,  the  coherence  length,  Dx>  ,  was 
the  distance  transverse  to  the  direction  of  propagation  where 
the  spatial  coherence  of  the  <E*E>  wave  dropped  to  e_1  . 

For  a  continuous  plot  of  the  atmospheric  MCF,  the 
e_1  distance  could  be  picked  off  the  curve  easily.  Due  to 
the  discrete  nature  of  the  simulation,  the  MCF  curve  was  not 
continuous,  but  a  series  of  points  equal  to  the  choice  of  the 
pixel  width  of  the  aperture  array.  Therefore  this  simulation 
used  a  linear  interpolation  routine  to  pick  out  values  for 
the  coherence  length  as  seen  in  the  results  of  the  following 
section . 
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IV.   RESULTS 

In  order  to  verify  that  this  numerical  simulation  produced 
correct  results,  we  had  to  find  ways  to  check  the  validity  of 
the  code.  In  the  code,  presented  in  Appendix  C,  there  are 
several  checkpoints  where  the  accuracy  of  the  simulation  could 
be  verified.  The  following  sections  discuss  these  checkpoints 
in  increasing  order  of  sophistication. 

A.   APERTURE  MCF 

An  initial  point  to  check  the  accuracy  of  the  simulation 
code  was  to  look  at  the  mutual  coherence  function  (MCF)  of 
the  aperture.  Since  the  code  offered  two  choices  of  aperture 
shape,  square  and  circular,  both  were  verified. 

Using  the  two-dimensional  FFT  routine,  the  diffraction 
patterns  for  100  mm  wide  square  and  circular  apertures  were 
calculated.  The  square  aperture  diffraction  pattern  in  Figure 
4.1a  displays  the  familiar  [(sin  x)/x]2  form  discussed  by 
Hecht  with  the  correct  maximum  value  and  node  spacing.  [Ref . 
5]  The  circular  aperture  diffraction  pattern  in  Figure  4.1b 
displays  the  Airy  pattern  derived  from  the  Ji  Bessel  function 
solution  carried  out  by  Hecht.   [Ref.  5] 

An  inverse  2-D  FFT  of  each  100  mm  wide  aperture  shape 
produced  the  aperture  MCF  (identical  to  aperture 
autocorrelation) .  The  square  aperture  autocorrelation  in 
Figure  4.2a   is   the  predicted  triangle  function.    The 
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Figure  4.1   Irradiance   Diffraction  Pattern  Resulting  From 
A  (a)  Square  Aperture  And  A  (b)  Circular  Aperture. 


33 


©      ©' 

uoi^Biajjooo^ny 


nor^Biejjooo^ny 


Figure  4.2   Autocorrelation  (Aperture  Mutual  Coherence 
Function)  Of  A  100  mm  Wide  (a)  Square  Aperture  And  A 
100  mm  Wide  (b)  Circular  Aperture 
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circular  aperture  autocorrelation  in  Figure  4.2b  is  as 
expected  for  a  "top  hat"  function. 

The  algorithm  gave  the  correct  MCF  results  which  verified 
the  simulation  up  through  the  aperture  MCF  point. 

B .  FILTERING 

The  most  crucial  portion  of  the  entire  simulation  code  was 
the  correct  frequency  filtering  algorithm.  A  succinct  way  of 
verifying  this  came  from  examining  the  result  of  taking  the 
logarithm  of  both  sides  of  Equation  3.3 

log{*f|  =  log{1.303(2n)Ml/X)2Cn*62f  -»W3) 
which  is 

log{*H  =  log{1.303(2n)3(l/X)2Cn252|  -  (11/3)  log{  f\  ,    (4.1) 
an  equation  for  a  straight  line  with  slope  of  -11/3. 

Figure  4.3  shows  a  plot  of  filtering  function  vs. 
frequency  f  with  a  value  10~13  for  Cn2  and  a  wavelength  of  600 
nm  (red  light)  using  the  Martin  and  Flatte  method  of  filtering 
discussed  earlier.  A  linear  fit  to  the  data  points  down  to 
the  Nyquist  frequency  had  a  slope  of  -3.88.  This  was  a 
deviation  from  the  predicted  value  of  -11/3  (-3.83)  by  1.3%. 
Therefore,  it  appears  that  the  filtering  is  being  accomplished 
correctly.  The  upturn  of  points  occurs  at  the  Nyquist  folding 
point  that  is  present  in  the  array  scheme  for  the  FFT  routine. 

C.  ATMOSPHERIC  MCF 

An  important  place  to  check  results  of  the  simulation  was 
at   the  calculation  of  the  atmospheric  MCF.   Following  the 
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Figure  4.3    Slice  Through  Phase  Screen  Function 
Filtered  By  The  Martin  And  Flatte  Method  Showing 

Its  Dependence  On  Spatial  Frequency. 
Slope  Of  -11/3  Comes  Directly  From  The  Kolmogorov 
Power  Law.   The  Upturn  At  Higher  Frequencies  Is  Due 
To  Nyquist  Folding  Present  In  The  Phase  Screen  Array  Code 
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procedure  for  calculating  the  atmospheric  MCF  outlined  in 
section  III,  first  the  diffraction  pattern  was  calculated. 
Figure  4.4a  shows  the  diffraction  pattern  for  a  100  mm  wide 
square  aperture  that  results  for  red  light  (A  =  0.6  urn)  with 
Cn2  =  10"16.  Notice  that  the  diffraction  is  only  slightly 
affected  from  what  is  shown  in  Figure  4.1a.  As  the  turbulence 
structure  parameter  Cn2  increased,  the  diffraction  pattern 
began  to  spread  out  more  until  it  no  longer  resembles  the 
[(sin  x)/x]2  form  it  had  without  turbulence,  as  seen  in 
Figures  4.4b,  4.4c,  and  4 . 4d  for  values  of  Cn2  =  10~18,  Cn2  = 
10~14,  and  Cn2  =  10~13  respectively.  Figure  4.4c  shows  that 
modest  amounts  of  turbulence  displace  the  image  centroid,  due 
to  low  frequency  tilting  of  the  electric  field.  Since  the 
diffraction  pattern  is  analogous  to  the  intensity  that  would 
be  seen  at  the  image  plane,  this  clearly  shows  the  blurring 
of  an  image  commonly  seen  looking  through  a  turbulent 
atmosphere  as  well  as  the  "dancing"  image  effect  at  lower 
levels  of  turbulence. 

The  atmospheric  MCF  was  calculated  for  different  values 
of  turbulence  by  the  method  described  in  section  III  using  the 
classical  method  of  filtering.  These  are  seen  in  Figures  4.5 
for  Cn2  values  of  0,  10"1*  ,  10~10,  10-1*  ,  and  10-13.  With  no 
turbulence,  the  atmospheric  MCF  is  a  value  of  1.0  everywhere, 
seen  as  a  straight  line.  With  increasing  turbulence,  the 
atmospheric  MCF  curves  fall  off  rapidly  to  zero. 
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(a) 


(b) 


(c) 


(d) 


Figure  4.4   Diffraction  Patterns  For  A  100  mm  Wide  Square 

Aperture  For  Four  Different  Values  Of  Turbulence: 

(a)  C2=10-»«  (b)  C2=10-i9  (c)  c„2=10-n  (d)  C„2=10-13. 

Displacement  Of  The  Beam  Centroid  Is  Obvious  In  (c)  Due 

To  Low  Frequency  Tilting  Of  The  Electric  Field. 
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Figure  4.5   Atmospheric  Mutual  Coherence  Function  Curves 
For  Five  Different  Levels  Of  Turbulence. 
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D.   COHERENCE  LENGTH 

The  e-1  point  of  atmospheric  MCF  curves  give  the  coherence 
length,  no.  In  the  simulation,  we  achieved  this  by  a  method 
of  linear  interpolation  as  discussed  previously.  One  should 
not  put  too  much  faith  in  a  statistical  sample  of  one, 
however.  For  this  reason,  a  modification  of  the  simulation 
was  made  to  run  several  samples  of  each  different  set  of  input 
parameters  so  that  a  better  statistical  representation  of 
results  occurred. 

Theoretical  coherence  lengths  were  calculated  from 
Equation  2.41  [Ref .  1]  for  three  different  values  of  Cn2  using 
R=1000  m  and  A=0.6  pm.  Several  runs  of  data  were  produced  to 
investigate  the  effect  of  different  sizes  of  both  the  electric 
field  array  and  the  aperture  array.  Calculated  results  are 
discussed  and  displayed  in  following  pages. 

Figure  4.6  shows  the  trend  for  coherence  length  values  as 
a  function  of  aperture  array  sample  points  using  the  classical 
method  of  filtering  for  Cn2  =  10-13.  Ten  samples  were  obtained 
for  each  data  point.  As  discussed  in  section  III,  aliasing 
makes  the  values  of  Do  too  high  if  there  are  too  few  sample 
points.  Using  aperture  arrays  no  larger  than  one-half  the 
width  of  the  main  array  was  done  to  avoid  tilt  problems.  The 
curves  asymptotically  approach  coherence  length  values  of  4.18 
mm  for  the  512  x  512  array  and  4.49  mm  for  the  256  x  256 
array.  These  represent  errors  of  39.3%  and  49.7% 
respectively. 
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Figure  4.6    Calculated  Coherence  Length  Values  For 

Cn2  =  10-13  Using  Different  Aperture  Array  Sizes 

And  Classical  Method  Of  Filtering. 

Error  Bars  Shown  Are  For  Case  Of  512  By  512  Array  Size, 

Representing  The  Standard  Deviation  For  Ten  Samples. 
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Figure  4.7  shows  coherence  length  trends  for  Cn2  = 
10_1«  also  using  the  classical  method  of  filtering.    Ten 
samples  were  obtained  for  each  data  point.  Results  were  15.78 
mm  for  the  512  x  512  array  and  18.43  mm  for  the  256  x  256 
array,  representing  errors  of  32.1%  and  54.2%  respectively. 

Figure  4.8  shows  coherence  length  trends  for  Cn2  = 
10-19  also  using  the  classical  method  of  filtering.    Ten 
samples  were  obtained  for  each  data  point.  Results  were  61.99 
mm  for  the  512  x  512  array  and  76.27  mm  for  the  256  x  256 
array,  representing  errors  of  23.3%  and  60.3%  respectively. 

Figure  4.9  shows  coherence  length  trends  for  Cn2  =  10-13 
using  the  Martin  and  Flatte  method  of  filtering.  Ten  samples 
were  obtained  for  each  data  point.  Results  were  4.07  mm  for 
the  512  x  512  array  and  4.46  mm  for  the  256  x  256  array, 
representing  errors  of  35.7%  and  48.7%  respectively. 

Figure  4.10  shows  coherence  length  trends  for  Cn2  =  10- l4 
also  using  the  Martin  and  Flatte  method  of  filtering.  Ten 
samples  were  obtained  for  each  data  point.  Results  were  14.94 
mm  for  the  512  x  512  array  and  20.46  mm  for  the  256  x  256 
array,  representing  errors  of  25.0%  and  71.2%  respectively. 

Figure  4.11  shows  coherence  length  trends  for  Cn2  =  10_1° 
also  using  the  Martin  and  Flatte  method  of  filtering.  Ten 
samples  were  obtained  for  each  data  point.  Results  were  60.74 
mm  for  the  512  x  512  array  and  74.04  mm  for  the  256  x  256 
array,  representing  errors  of  27.7%  and  55.6%  respectively. 
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Figure  4.7    Calculated  Coherence  Length  Values  For 

Cn2  =  10-14  Using  Different  Aperture  Array  Sizes 

And  Classical  Method  Of  Filtering. 

Error  Bars  Shown  Are  For  Case  Of  512  By  512  Array  Size, 

Representing  The  Standard  Deviation  For  Ten  Samples. 
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Figure  4.8    Calculated  Coherence  Length  Values  For 

Cn2  =  10-1B  Using  Different  Aperture  Array  Sizes 

And  Classical  Method  Of  Filtering. 

Error  Bars  Shown  Are  For  Case  Of  512  By  512  Array  Size, 

Representing  The  Standard  Deviation  For  Ten  Samples. 
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Figure  4.9    Calculated  Coherence  Length  Values  For 

Cn2  =  10-13  Using  Different  Aperture  Array  Sizes 

And  Martin  &  Flatte  Method  Of  Filtering. 

Error  Bars  Shown  Are  For  Case  Of  512  By  512  Array  Size, 

Representing  The  Standard  Deviation  For  Ten  Samples. 
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Figure  4.10    Calculated  Coherence  Length  Values  For 
Cn   -  10-»«  Using  Different  Aperture  Array  Sizes 
And  Martin  &  Flatte  Method  Of  Filtering 
Error  Bars  Shown  Are  For  Case  Of  512  By  512  Array  Size 
Representing  The  Standard  Deviation  For  Ten  Samples 
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Figure  4.11    Calculated  Coherence  Length  Values  For 

Cn2  =  10"1B  Using  Different  Aperture  Array  Sizes 

And  Martin  &  Flatte  Method  Of  Filtering. 

Error  Bars  Shown  Are  For  Case  Of  512  By  512  Array  Size, 

Representing  The  Standard  Deviation  For  Ten  Samples. 
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The  most  accurate  results  were  obtained  when  the  aperture 
array  was  one-fourth  the  size  of  the  electric  field  array,  as 
seen  in  Figures  4.6  through  4.11.  For  this  reason  all 
subsequent  coherence  lengths  were  calculated  using  a  512  by 
512  array  with  a  128  by  128  aperture  array. 

Figure  4.12  shows  coherence  length  values  calculated  with 
different  values  of  turbulence  using  the  classical  method  of 
filtering.  Ten  samples  were  obtained  for  each  calculated 
value.  Although  the  calculated  coherence  length  values 
deviated  from  theoretical  values  by  approximately  30%,  a 
linear  fit  to  the  data  shows  that  the  slopes  were  correct, 
implying  a  constant  error  in  the  algorithm  due  to 
underestimating  turbulence  effects. 

Adding  a  log-normal,  random-amplitude  screen  to  the  GAUSS 
subroutine  was  attempted  to  improve  the  simulation  results  for 
coherence  lengths.  Having  both  phase  and  amplitude  screens 
present  is  equivalent  to  the  complete  Rytov  approximation, 
Equation  2.31.  Eight  samples  were  obtained  for  each 
calculated  value.  Figure  4.13  dramatically  shows  that  this 
simple  addition  produced  calculated  coherence  lengths  that 
were  within  3%  of  theoretical  values,  from  a  linear  fit  to 
the  data. 

Table  4.1  presents  coherence  lengths  calculated  by  the 
simulation  using  both  random  phase  and  amplitude  screens  in 
the  complete  Rytov  approximation.   All  data  points  were 
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Figure  4.12    Calculated  Coherence  Length  Values  As  A 
Function  Of  Cn2  With  Only  Random  Phase 
Screen  In  The  Simulation. 
Calculated  Values  Are  About  30%  High  For  All  Levels  Of 
Turbulence,  From  Linear  Fit  To  The  Data. 
Error  Bars  Represent  Standard  Deviation  For  Ten  Samples 
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Figure  4.13   Calculated  Coherence  Length  Values  As  A 
Function  of  Cn2  With  Both  Random  Phase  And  Amplitude 

Screens  In  The  Simulation. 
Calculated  Values  Are  About  3%  High  For  All  Levels  Of 
Turbulence,  From  Linear  Fit  To  The  Data. 
Error  Bars  Represent  Standard  Deviation  For  Eight  Samples 
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obtained  using  a  512  by  512  array  and  128  by  128  aperture 
array. 

TABLE  4.1   COHERENCE  LENGTHS 

This  table  shows  theoretical  values  of  coherence  lengths 
calculated  from  Equation  2.39  for  eight  different  values  of 

Cn2  . 


Cn2 

theoretical  Do     (mm) 

calculated  po    (mm) 

1    X    lO"19 

47.57 

50        *    30 

2   x   10-19 

31.39 

30        *    20 

5   x    10-13 

18.11 

14        *      2 

1    x   10*14 

11.95 

9*3 

2   x    10"14 

7.88 

7*2 

5   x   10-n 

4.55 

4*1 

1    x   10-13 

3.00 

3*1 

2   x   10-13 

1.98 

2.0    *      0.4 

E.   INTENSITY  VARIANCE  SATURATION 

The  last  checkpoint  available  for  verifying  proper 
operation  of  the  simulation  was  whether  the  code  led  to  the 
saturation  of  intensity  variance  for  increasing  turbulence. 
This  should  happen  because  the  Rytov  assumption  of  an 
exponential  random  term,  exp(^),  has  a  magnitude  bounded  by 
plus  or  minus  one.  Figure  4.14  shows  that  saturation  does 
occur  around  Cn2  =  10-14  at  a  value  of  unity,  as  expected.  The 
normalized  variance  >  1  present  at  this  level  of  turbulence 
is  due  to  a  "dancing"  beam  centroid  caused  by  low  frequency 
tilt  of  the  electric  field  as  discussed  previously.  A  decline 
to  unity  normalized  variance  appears  in  the  calculations  for 
higher  turbulence  values.   [Ref.  12] 
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Figure  4.14    Normalized  Intensity  Variance  As 
A  Function  Of  Cn2  . 
Saturation  Occurs  Around  Cn2  =  10~14  . 
Large  Variances  Near  Saturation  Are  Due  To  Beam 
Centroid  Displacement  Cause  By  Low  Frequency  Tilt. 
Error  Bars  Represent  Standard  Deviation  For  Ten  Samples 
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V.   CONCLUSIONS 

The  program  coding  discussed  in  this  thesis  simulated  the 
propagation  of  optical  waves  through  a  random  media  using  two- 
dimensional  fast  Fourier  transform  (FFT)  techniques.  The 
coding  used  Fraunhofer  propagation  throughout.  An  extension 
to  include  multi-step  Fresnel  propagation  is  straightforward. 

Aliasing  problems  in  the  FFT  were  avoided  by  taking  sample 
points  closer  together  than  the  theoretical  coherence  length 
for  a  given  turbulence  structure  parameter,  Cn2  .  Tilt 
problems  in  the  FFT  were  suppressed  by  choosing  the  aperture 
array  no  larger  than  one-half  the  size  of  the  FFT  array,  with 
one-fourth  size  giving  the  best  results. 

Accuracy  of  the  simulation  was  verified  at  various 
checkpoints  within  the  coding.  Aperture  diffraction  patterns 
and  autocorrelations  were  compared  to  expected  results. 
Turbulence  effected  diffraction  patterns  and  atmospheric 
mutual  coherence  functions  were  presented  to  show  trends  for 
increasing  turbulence.  Calculated  coherence  length  values 
were  approximately  30%  larger  than  theoretical  coherence 
values  when  using  a  512  by  512  FFT  array  with  a  128  by  128 
aperture  array.  These  results  were  obtained  with  only  a 
Gaussian  distributed  random  phase  screen  in  the  simulation. 

The  addition  of  a  Gaussian  distributed  random  amplitude 
screen   (for  the  complete  Rytov  approximation)    to   the 
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simulation  brought  calculated  coherence  values  to  within  3% 
of  theoretical  values.  This  showed  that  a  numerical  error  was 
not  present  in  the  coding,  but  that  turbulence  effects  were 
underestimated  without  the  random  amplitude  term  included. 
A  single  step,  Fraunhofer  algorithm  must  include  the  amplitude 
term  to  simulate  turbulence  effects  on  an  optical  wave 
correctly.  Unity  saturation  of  the  normalized  intensity 
variance  occurred  for  Cn2  values  of  approximately  10-14. 

Coherence  length  values  used  two  types  of  filtering.  The 
classical  method  and  the  Martin  and  Flatte  method  gave 
comparable  results.  The  computational  efficiency  of  the 
Martin  and  Flatte  is  significant. 

Future  endeavors  with  this  simulation  should  include  the 
incorporation  of  Fresnel  propagation  into  the  coding.  Also, 
further  studies  with  the  complete  Rytov  approximation  in  the 
simulation  code  should  be  pursued  to  obtain  better  statistical 
results . 
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APPENDIX  A.   RANDOM  NUMBER  GENERATOR  SUBROUTINE 


subroutine  RAN(iran,r) 
c 

c   This  is  the  random  number  generator  algorithm  from  Dr. 
c   Harrison's  notes,  (Ref.  7). 
c 

c   Variables: 
c     iran  —  input  seed  value  (5  digit  integer) 

c     r returned  random  number  uniformly  distributed  from  0  to  1 

c 

iran=iran*99947 

r=0.5  +  real(iran)*2.328306e-10 
c 

return 

end 
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APPENDIX  B.   GAUSSIAN  DISTRIBUTION  SUBROUTINE 


subroutine  GAUSS (narray,iran) 
c 

c   This  subroutine  changes  uniformly  distributed  random  numbers, 
c   provided  by  subroutine  RAN,  to  Gaussian  distributed  random 
c   numbers  with  unity  variance.  This  subroutine  is  found  in 
c   Knuth,  [Ref.  8]. 
c 
c   Variables: 

c     phaser real  phase  screen  2-D  array 

c     phasei imaginary  phase  screen  2-D  array 

c     narray dimension  of  the  2-D  arrays 

c     iran input  seed  value  for  RAN  (5  digit  integer) 

c     r uniformly  distributed  random  numbers 

c     xl  &  x2 Gaussian  distributed  random  numbers 

c 

common  /blk2/phaser (512,512) , phasei (512,512) 
c 

do  10  i=l, narray 
do  10  j=l, narray 
20     call  RAN(iran,r) 
vl=2.*r-l. 
call  RAN(iran,r) 
v2=2.*r-l. 
s=vl*vl  +  v2*v2 
if  (s.ge.1.0)  goto  20 
scale=sqrt (-2.*alog(s)/s) 
xl=vl*scale 
x2=v2*scale 
phaser (i,j)=xl 
phasei (i, j)=x2 
10  continue 
c 

return 
end 
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APPENDIX  C.   SIMULATION  CODE 


c  This  code  simulates  the  propagation  of  a  monochromatic  optical 

c  wave  in  a  turbulent  atmosphere.  Tiro-dimensional  FFT  routines 

c  are  used  extensively  for  calculations.  Only  Fraunhofer 

c  propagation  is  present  in  the  coding.  Filtering  is  accomplished 

c  by  a  choice  of  two  methods:  (1)  classical  and  (2)  the  Martin 

c  and  Flatte  method. 

c 

c  Subroutines: 

c     INIT Sets  electric  field  arrays  to  zero.  A  value  of  0 

c  in  the  call  initializes  real  and  imaginary  arrays, 

c  A  value  of  1  initializes  only  the  imaginary  array, 

c  SQUARE. ..Establishes  the  planar  electric  field  for  the  case 

c  of  a  square  aperture. 

c  CIRCLE  —  Establishes  the  planar  electric  field  for  the  case 

c  of  a  circular  aperture. 

c     PLOT Gives  a  screen  plot  of  the  array  chosen  in  the  call 

c  using  EGA  graphics.  Different  colors  are  chosen  to 

c  represent  different  orders  of  magnitude  values.  Most 

c  calls  within  this  subroutine  are  NDP  Fortran-386 

c  specific. 

c     XFORM Takes  the  2-D  FFT  of  the  two  arrays  given  in  the  call. 

c  Returns  values  in  the  same  arrays.  Makes  calls  to 

c  subroutine  FFT  (1-D  FFT  routine)  several  times.  A 

c  value  of  -1  in  the  call  is  for  a  direct  transform, 

c  A  value  of  +1  in  the  call  is  for  an  inverse  transform. 

c     FFT 1-D  FFT  routine  supplied  by  Dr.  Don  Walters.  This 

c  routine  is  used  by  subroutine  XFORM  multiple  times 

c  to  accomplish  the  2-D  transform, 

c  This  is  Dr.  Walters  FFT 

c     MAG Calculates  the  electric  field  magnitude  and  the 

c  intensity  field  from  the  real  and  imaginary  electric 

c  field  arrays. 

c  MCFPLOT. .Calculates  the  Mutual  Coherence  Function  (MCF)  from 

c  the  intensity  field  and  then  plots  it  on  the  screen 

c  using  EGA  graphics.  Many  calls  within  this 

c  subroutine  are  NDP  Fortran-386  specific. 

c     GAUSS Creates  the  phase  screens  from  Gaussian  distributed 

c  random  numbers  with  unit  variance.  Makes  calls 

c  to  subroutine  RAN  for  uniformly  distributed  random 

c  numbers  then  transforms  them  to  a  Gaussian  distribution, 

c  Elements  of  this  subroutine  provided  by  Knuth.  [Ref.  8] 

c     RAN Generates  uniformly  distributed  random  numbers.  This 

c  algorithm  from  Harrison  [Ref.  7] 
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c     FILTER. . .Filters  the  phase  screens  according  to  the  Kolmogorov 
c  power  law.  A  shuffling  of  array  values  is  necessary 

c  for  proper  filtering. 

c     MESH Meshes  the  random  phase  screen  with  the  electric 

c  field. 

c 

c   Variables: 

c     fieldr. . .Real  2-D  electric  field  array. 

c     fieldi. . .Imaginary  2-D  electric  field  array. 

c     fieldro. .Real  2-D  planar  electric  field  at  aperture. 

c     phaser Real  2-D  phase  screen  array. 

c     phasei Imaginary  2-D  phase  screen  array. 

c     field*... 2-D  array  that  represents  the  electric  field  magnitude. 
c     fieldm2..2-D  array  that  represents  the  intensity  field, 
c     nar ray. . .Dimension  of  electric  field  array 

c     nfield Dimension  of  aperture  array 

c     iran Input  seed  for  random  number  generator  (5  digit  integer), 

c     cn2 Refractive  index  structure  parameter. 

c     wvl Wavelength  of  monochromatic  electromagnetic  wave. 

c     delz Physical  distance  between  aperture  plane  and  image  plane, 

c     width Physical  width  of  aperture  plane. 

c     delx Physical  distance  between  sample  points  in  aperture 

array. 

c     rhonot Coherence  length  calculated  by  interpolation  between 

c  points  of  MCF  curve. 

c     re Real  1-D  array  used  by  subroutine  FFT. 

c     rim Imaginary  1-D  array  used  by  subroutine  FFT. 

c     f phaser. .Real  2-D  array  used  in  subroutine  FILTER  that  represents 

c  the  phase  screen  in  a  more  convenient  form  for 

c  filtering  and  viewing. 

c     f phasei. .Imaginary  2-D  array  used  in  subroutine  FILTER  that 

c  represents  the  phase  screen  in  a  more  convenient  form 

c  for  filtering  and  viewing. 

c 

c 

common  /blkl/fieldr (512, 512) ,f ieldro(512, 512) ,fieldi (512,512) 

common  /blk2/phaser (512,512) , phasei (512,512) 

common  /blk3/fieldm(512,512) ,f ieldm2(512,512) 
c 

c   Initialize  arrays 
c 

1  call  INIT(narray,0) 
c 

c   Input  section 
c 

write (*,*)  'Enter  dimension  of  array  that  required.' 

read(*,*)  narray 

write(*,*)  'Enter  dimension  of  planar  electric  field.' 

write(*,*)  '(Pixel  width  of  the  aperture)' 

read(*,*)  nfield 

2  write (*,*)  "Choose  aperture  shape:   1)  SQUARE' 
write(*,M  '  2)  CIRCLE' 
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read(*,*)  ichoice 

if  (ichoice  .It.  1  .or.  ichoice  .gt.  2)  then 

write(*,*)  'Try  again!' 

goto  2 
endif 

write(*,*)  'Enter  the  seed  for  the  random  number  generator.' 
write(*,*)  '(Value  mist  be  a  five  digit  integer)' 
read(*,*)  iran 

write (*,*)  'Enter  the  value  of  Cn  squared.' 
read(*,*)  cn2 

write(*,*)  'Enter  the  wavelength  of  light.' 
read(*,*)  wvl 

write(*,*)  'Enter  the  distance  from  aperture  to  observer.' 
read(*,*)  delz 
delz=1000. 

write(*,*)  'Enter  the  width  of  the  aperture  in  meters.' 
read(*,*)  width 


delx=width/real (nf ield) 
c 

c   Load  arrays  with  desired  planar  electric  field 
c 

if  (ichoice  .eq.  1)  call  SQUARE (narray,nf ield) 

if  (ichoice  .eq.  2)  call  CIRCLE (narray,nf ield) 
c 

c   Plot  planar  electric  field 
c 

call  PLOT(fieldr,narray,l) 
c 

c   Take  Fourier  transform  of  the  field 
c 

call  XFORM(f ieldr,f ieldi,narray,delx,-l.) 
c 

c   Calculate  magnitude  of  transformed  field  then  plot  it 
c 

call  HAG(narray) 
c 

call  PLOT(f ieldm,narray,l) 
c 

c   Set  imaginary  portion  of  field  to  zero 
c 

call  INIT(narray,l) 
c 

c   Take  Fourier  transform  of  field  intensity 
c 

call  XF0RM(fieldm2,fieldi,narray,delx,+l.) 
c 

c   Calculate  and  plot  MCF  of  aperture 
c 

call  MCFPLOT(narray,delx,nfield,0) 
c 
c   Load  phase  screen  arrays  with  Gaussian  random  numbers 
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call  GAUSS (narray,iran) 
c 

c   Transform  phase  screen  to  frequency  space 
c 

call  XFORM (phaser , phasei , narray , sqrt (real (narray ) ) , -1 . ) 
c     call  PLOKphaser, narray,  1) 
c     call  PLOTtphasei, narray, 1) 
c 

c   Filter  phase  screen  using  Kolmogorov  spectrum  idea 
c 

call  FILTER (narray , cn2 , wvl , delx , delz) 
c     call  PLOT(phaser, narray, 1) 
c     call  PLOT(phasei, narray, 1) 
c 

c   Transform  phase  screen  back  to  real  space 
c 

call  XFORM (phaser, phasei, narray, sqrt (real (narray )*delx) ,+1.) 
c     call  PLOT (phaser, narray, 1) 
c     call  PLOT (phasei, narray, 1) 
c 

c   Mesh  phase  screen  Kith  planar  electric  field 
c 

call  INIT(narray,l) 

call  MESH (narray) 
c 

call  PLOT(fieldr, narray, 1) 
c 

c   Take  fourier  transform  of  this  electric  field 
c 

call  XFORM ( f ieldr ,f ieldi, narray, delx, -1.) 
c 

call  MAG (narray) 
c 

call  PLOT(fieldm, narray,  1) 
c 

c   Reset  imaginary  portion  of  field  to  zero 
c 

call  INIT(narray,l) 
c 

c   Take  Fourier  transform  to  get  diffraction  pattern  w/  turbulence, 
c 

call  XFORM (fieldm2,f ieldi, narray, delx, +1.) 
c 

c   Calculate  atmospheric  MCF  and  plot  it. 
c 

call  MCFPLOT (narray, delx, nfield,l) 
c 

goto  1 

stop 

end 
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subroutine  INIT(narray,k) 


common  /blkl/fieldr (512,512) ,fieldro(512, 512) ,fieldi(512, 512) 
c     do  10  i=l,narray 

do  10  j=l,narray 
if  (k  .eq.  0)  then 
fieldr(i,j)=0.0 
fieldro(i, j)=0.0 
endif 

£ieldi(i, j)=0.0 
10  continue 
c 

return 
end 


subroutine  MAG(narray) 

common  /blkl/fieldr (512,512) ,fieldro(512, 512) ,fieldi  (512,512) 
common  /blk3/f ieldm(512,512) ,f ieldm2(512,512) 

do  10  i=l,narray 
do  10  j=l,narray 

fieldm2(i,j)=fieldr(i,j)**2  +  f ieldi(i,j)**2 
f ieldm(i, j)=sqrt (f ieldm2(i, j) ) 
10  continue 

return 
end 

subroutine  SQUARE (narray,nf ield) 

common  /blkl/fieldr (512,512) ,fieldro(512, 512) ,fieldi (512,512) 

do  10  i=narray/2+l-nfield/2,narray/2+nfield/2 
do  10  j=narray/2+l-nfield/2,narray/2+nfield/2 
fieldr(i,j)=1.0 
fieldro(i, j)=1.0 
10  continue 

return 
end 

subroutine  CIRCLE(narray,nf ield) 

common  /blkl/fieldr (512, 512) , fieldro(512, 512) ,fieldi(512, 512) 

narray2=narray/2 

nfield2=nfield/2 

do  10  i=narray2+l-nf ield2,narray2+nf ield2 
do  20  j=narray2+l-nf ield2,narray2+nf ield2 
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radius2= (real (i) -real (narray2) ) **2  + 
*  (real (j) -real (narray2))**2 

if  (radius2  .It.  reaHnf  ield2)**2)  then 
fieldr(i,j)=1.0 
fieldro(i,j)=1.0 
endif 
20   continue 
10  continue 
c 

return 
end 

c 

subroutine  PLOT(field,narray,k) 
c 

dimension  f ield(512,512) 
dimension  ndex(20) 

data  ndex/4, 4 ,12, 12, 14, 14, 10, 10, 2, 2, 3, 3, 9, 9, 1,1 ,8, 8,0,0/ 
c 

fmax=0.0 

do  100  i=l,narray 
do  100  j=l,narray 
x=field(i, j) 
if  (x.gt.fnax)  then 

faax=x 
endif 
100  continue 
c 

call  set_video_node(16) 
call  ega_set_mode_4 
c 

do  110  i=l,narray/k 
do  110  j=l,narray/k 
ix=i 

iy=j 

xmax=alogl0 (field (i , j ) It  max) 
index=8*abs (xmax) +1 
if  (index. ge. 21)  then 

icolor=0 
else 

icolor=ndex (index) 
endif 

call  ega_put_pixel (ix,iy,icolor) 
110  continue 
c 

call  pause 

call  set_video_mode(3) 
c 

return 
end 

c 

subroutine  MCFPLOT(narray,delx,nf ield,iaprture) 
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common  /blk3/f ieldm(512,512) ,fieldm2(512,512) 

dimension  apmcf (513) ,ymcf (513) 

logical  flag 

efold=0. 36787944 

rhonot=0. 

£lag=.true. 

do  10  i=l,513 

if  (iaprture  .eq.  0)  then 

apmcf (i)=0.0 
endif 

ymcf (i)=0.0 
10  continue 

fmcf=0.0 

do  20  i=l,narray 
do  20  j=l,narray 
xmcf=fieldm2(i, j) 
if  (xmcf .gt.fmcf )  then 

fmcf=xmcf 
endif 
20  continue 

do  30  i=l,nfield 

if  (iaprture  .eq.  0)  then 
apmcf (i) =f ieldm2 (1 , i) /f mcf 
ymcf (i)=apmcf (i) 
else 

ymcf (i)=f ieldm2(l,i)/ (apacf (i)*f mcf) 
endif 
30  continue 

do  40  i=l,nfield 

if  (iaprture. eq.l)  then 
if  (flag)  then 

if  (ymcf (i).lt.efold)  then 
rhonot=real (i-1) *delx 
*  +  (efold-ymcf (i))*delx/ (ymcf (i-1) -ymcf (i)) 

flag=. false, 
endif 
endif 
endif 
40  continue 

call  set_video_mode(16) 
call  ega_set_mode_4 
call  locate(20,l) 
if  (iaprture  .eq.  0)  then 

call  write_string( 'Aperture  MCF  vs.  Length') 
else 

call  write_string( 'Atmospheric  MCF  vs.  Length') 
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endif 

call  located, 1) 
call  *rite_string('1.0' ) 
call  locate(25,24) 
call  write_string (' length ' ) 
c     call  locate(37,24) 
c     call  write_string(nfield) 

call  ega_draw_line(40,20,40,320,15) 
call  ega_drav_line(40,320,600,320,15) 
c 

ix=40 
iy=20 

do  50  i=l,nfield 
ix2=40+560*i/nfield 
iy2=20+int (300.* (l.-ymcf (i+1) ) ) 
call  ega_draw_line ( ix, iy ,  ix2 , iy2 , 15 ) 
ix=ix2 
iy=iy2 
50  continue 
c 

call  pause 

call  set_video_mode(3) 

write(*,*)  'The  coherence  length  is  ',rhonot,'  meters. ' 

return 
end 

subroutine  XF0RM(f ieldr, fieldi ,nar ray, fftnorm, sign) 

common  /blk4/re(512) ,rim(512) 

dimension  fieldr (512,512) ,f ieldi(512,512) 

data  re/512*0./,rim/512*0./ 

m=int (alog (real (narray) ) /alog (2 . ) ) 

do  30  i=l, narray 
do  40  j=l, narray 
re(j)=fieldr(i,j) 
rim(j)=fieldi(i, j) 
40   continue 

call  FFT(m,f f tnorm,sign) 
do  50  j=l, narray 
fieldr(i, j)=re(j) 
fieldi(i, j)=rim(j) 
50   continue 
30  continue 

do  60  j=l, narray 
do  70  i=l, narray 
re(i)=fieldr(i, j) 
rim(i)=f ieldi(i, j) 
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70 

continue 

call  FFT(m,fftnorm,sign) 

do  80  i=l,narray 

fieldrd,  j)=re(i) 

fieldi(i, j)=rim(i) 

80 

continue 

60 

continue 

return 

end 

subroutine  FFT(m,ff tnorm,sign) 

common  /blk4/re(512) ,rim(512) 

pi=3.141592653589792*sign 

n=2**m 

nl=n-l 

3=1 

do  200  i=l,nl 
if  (i.lt.j)  then 
t=re(j) 
re(j)=re(i) 
re(i)=t 
t=rim(j) 
rim(j)=rim(i) 
rim(i)=t 
endif 
k=n/2 

do  201  while  (k.lt.j) 
j-j-k 
k=k/2 
201   continue 

j=3+k 
200  continue 
le=l 

do  202  1=1, m 
lel=le 
le=le+le 
ure=l. 
uim=0. 
ang=pi/lel 
wre=cos(ang) 
wim=sin(ang) 
do  203  j=l,lel 
do  204  i=j,n,le 
ip=i+lel 

tre=re (ip) *ure-rim (ip) *uim 
tim=re (ip) *uim+rim (ip) *ure 
re(ip)=re(i)-tre 
rim(ip)=rim(i)-tim 
re(i)=re(i)+tre 
rim(i)=rim(i)+tim 
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204  continue 
t=ure*wre-uim*wim 
uim=ure*wim+uim*wre 
ure=t 

203   continue 
202  continue 

if  (sign. gt. 0.0)  then 

pts=1.0/real(n*f f tnorm) 
else 

pts=ff tnorm 
endif 
do  205  i=l,n 

re(i)=re(i)*pts 

rim(i)=rim(i)*pts 

205  continue 

return 
end 


subroutine  GAUSS (narray , iran) 

common  /blk2/phaser (512,512) ,phasei(512,512) 

do  10  i=l, narray 
do  10  j=l,narray 
20     call  RAN(iran,r) 

vl=2.*r-l. 

call  RAN(iran,r) 

v2=2.*r-l. 

s=vl*vl  +  v2*v2 

if  (s.ge.1.0)  goto  20 

scale=sqrt (-2.*alog(s) /s) 

xl=vl*scale 

x2=v2*scale 

Phaser (i, j)=xl 

phasei (i, j)=x2 
10  continue 

return 
end 

subroutine  RAN(iran,r) 

kran=99947 

iran=iran*kran 

r=0.5  +  real(iran)*2.328306e-10 

return 
end 

subroutine  FILTER (narray , cn2 , wvl , delx, delz) 
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common  /blk2/phaser (512,512) ,phasei(512,512) 

dimension  fphaser (512,512) ,fphasei(512,512) 

dimension  aaa(512,512) 

pi=3. 141592653589792 

tpi=2.*pi 

narray2=narray/2 

npivot=narray2+l 

power=-ll./6. 

f actor=sqrt (1 . 303* (tpi**3) *cn2*delz) /wvl 

factor2=(tpi/real(delx*narray)  )**power 

do  10  i=l,narray 
do  10  j=l,narray 
fphaser (i,j)=0. 
fphasei(i, j)=0. 
10  continue 

do  20  i=l,narray 
do  20  j=l,narray 

if  (i.le.npivot)  then 
if  (j .le.npivot)  then 

fphaser (i-l+narray2 , j-l+narray2) =phaser (i , j ) 
fphasei (i-l+narray2, j-l+narray2)=phasei (i, j) 
else 

fphaser (i-l+narray2, j-l-narray2)=phaser (i, j) 
fphasei (i-l+narray2, j-l-narray2)=phasei (i, j) 
endif 
else 

if  (j. le.npivot)  then 

fphaser (i-l-narray2, j-l+narray2)=phaser (i, j) 
fphasei (i-l-narray2, j-l+narray2)=phasei (i, j) 
else 

fphaser (i-l-narray2, j-l-narray2)=phaser (i, j) 
fphasei (i-l-narray2, j-l-narray2)=phasei(i, j) 
endif 
endif 
20  continue 

do  30  i=l,narray 
do  30  j=l,narray 

freq=sqrt(real(i-narray2)**2  +  real(j-narray2)**2) 
if  (i.eq.narray2.and. j .eq.narray2)  then 
fphaser (i,j)=0. 
fphasei(i, j)=0. 
else 

fphaser (i, j )=f phaser (i, j)*f actor *f actor 2* (f req) **power 
f phasei(i, j)=fphasei(i, j) *f actor *f actor 2* (f req)**power 
endif 
30  continue 
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c 

do  40  i=l,narray 
do  40  j=l,narray 

if  (i.lt .narray2)  then 
if  (j .lt.narray2)  then 

Phaser (i+l+narray2 , j+l+narray2) =f phaser (i , j ) 
phasei (i+l+narray2, j+l+narray2)=f phasei (i,j) 
else 

Phaser (i+l+narray2 , j+l-narray2) =f phaser (i , j ) 
phasei (i+l+narray2, j+l-narray2)=f phasei (i, j) 
endif 
else 

if  (j.lt.narray2)  then 

phaser ( i+l-narray2 , j+l+narray2) =f phaser ( i , j ) 
phasei ( i+l-narray2 , j+l+narray2) =f phasei (i , j ) 
else 

phaser (i+l-narray2, j+l-narray2)=f phaser (i,j) 
phasei (i+l-narray2, j+l-narray2) =f phasei (i , j ) 
endif 
endif 
40  continue 
c 

return 
end 

c 

subroutine  MESH(narray) 
c 

common  /blkl/fieldr (512,512) ,fieldro(512, 512) ,f ieldi(512,512) 
common  /blk2/phaser (512,512) , phasei (512,512) 
c 

do  10  i=l,narray 
do  10  j=l,narray 

ercosphi=f ieldro (i , j ) *cos (phaser (i , j ) ) *exp (phasei (i , j ) ) 
eicosphi=f ieldi (i, j) *cos (phaser (i, j) ) *exp (phasei (i, j) ) 
ersinphi=f ieldro ( i, j)*sin (phaser (i, j) )*exp (phasei (i, j) ) 
eisinphi=f ieldi (i, j) *sin (phaser (i, j) ) *exp (phasei (i, j) ) 
f ieldr (i , j ) =ercosphi-eisinphi 
f ieldi (i, j)=ersinphi+eicosphi 
10  continue 
c 

return 
end 
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